%% Generate random GP functions (one-dimensional)
% choose covar function
covfunc = 11;
switch covfunc
    case 1; 
        k = @(x,y) 1*x'*y;                      % Linear
    case 2; 
        k = @(x,y) 1*min(x,y);                  % Brownian motion
    case 3;
        k = @(x,y) exp(-100*(x-y)'*(x-y));      % squared exp
    case 4;
        k = @(x,y) exp(-1*sqrt((x-y)'*(x-y)));  % Ornstein-Uhlenbeck
    case 5;
        k = @(x,y) exp((-2*sin(1*pi*(x-y))^2));   % Periodic
    case 6;
        k = @(x,y) exp(-100*min(abs(x-y),abs(x+y))^2);  % A Symmetric GP
    case 7;
        k = @(x,y) exp(-2*sin(1*pi*(x-y))^2)+.2*exp(-2*sin(0.2*pi*(x-y))^2);   % Sum of Periodics (it works)
    case 8;
        k = @(x,y) exp(-1*sin(1*pi*(x-y) + 0.1*sin(0.4*pi*(x-y)))^2);
    case 9;
        k = @(x,y) exp(-1*(sin(1*pi*(x-y) + 0.5*(cos(0.4*pi*x)-cos(0.4*pi*y))))^2);
    case 10;
        k = @(x,y) exp((-8*cos(1*pi*(x-y))*(sin(0.4*pi*x+(pi/4)))^2* ...
            (sin(0.4*pi*y+(pi/4)))^2 + 4*(sin(0.4*pi*x+(pi/4)))^4 + 4*(sin(0.4*pi*y+(pi/4)))^4));
    case 11;
        k = @(x,y) exp(-2*(sin(.6*pi*(x-y))^2)/1^2)* exp(-1*(x-y)'*(x-y)/(5^2)); % Periodic + SE
end

% Choose points at which to sample
x = 0:0.05:10;
n = length(x);

% construct the covariance matrix
C = zeros(n,n);
for i = 1 : n
    for j = 1 : n
        C(i,j) = k(x(i),x(j));
    end 
end

% Sample from a Gaussian process at these points
u = randn(n,1);                 % sample u ~ N(0,I)
[A,S,B] = svd(C);               % factor C = ASB*
z = A*sqrt(S)*u;                % z = A S^.5 u ~ N(0,C) 

% Display function
figure; hold on; clf;
plot(x,z,'.-');
% axis([0,1,-2,2]);
